;+
; NAME:
;       IMPORT_DENSEGAS_SPECTRA
;
; PURPOSE:
;       Import the various dense gas tracer spectra into a single
;       velocity structure for analysis.
;
; CATEGORY:
;       distance-omnibus utility
;
; CALLING SEQUENCE:
;       IMPORT_DENSEGAS_SPECTRA [,/PLOT][,/HCOP][,/N2HP][,/NH3]
;                               [,/CS21][,ORDER=order][,/CREATE]
;
; INPUTS:
;       NONE  (reads in velocity IDL structure)
;
; OPTIONAL INPUTS:
;       ORDER  -- Order of polynomial to subtract from HHT HCO+ / N2H+
;                 spectra [Default: 15]
;
; KEYWORD PARAMETERS:
;       PLOT         -- Display plots of the input spectra as they are
;                       being processed.  Only displays plots of the
;                       line data selected with the [Line Names]
;                       keywords.
;       CREATE       -- Force creation of the IDL velocity structure.
;       [Line Names] -- Used with the /PLOT keyword, displays plots
;                       for this spectral line.
;
; OUTPUTS:
;       Places all the spectra into the velocity IDL structure
;
; OPTIONAL OUTPUTS:
;       NONE
;
; MODIFICATION HISTORY:
;
;       Created:  09/14/10, TPEB -- Initial Version
;       Modified: 10/19/10, TPEB -- Making polyfit order an adjustable
;                                   parameter
;       Modified: 07/08/11, TPEB -- Conformal change of PRO name
;                                   v_from_hdr()
;       Modified: 07/12/11, TPEB -- Offloaded baselining of spectra to
;                                   new routine BASELINE_SPEC.pro
;       Modified: 07/14/11, TPEB -- Cleaned up the memory hog by
;                                   undefining the big variables at
;                                   the end of the routine.  Also,
;                                   moved the GRS 13CO import to here
;                                   from the routine
;                                   create_velocity_structure.pro.
;       Modified: 08/30/11, TPEB -- Added spectrum baselining (via
;                                   BASELINE_SPEC.pro) to the code for
;                                   the CS transitions (2-1 and 5-4).
;       Modified: 10/17/11, TPEB -- Added spectrum filename
;                                   information to the velocity
;                                   structure for backwarks checking.
;                                   Also added /CREATE keyword to
;                                   force creation of new velocity
;                                   structure.
;       Modified: 10/18/11, TPEB -- Code cleanup.
;       Modified: 10/21/11, TPEB -- More code cleanup, placing reused
;                                   code sections into their own
;                                   functions for ease of updating
;                                   specifics.
;       Modified: 11/08/11, TPEB -- Removing use of the CSO CS(5-4)
;                                   sprectra from use in the
;                                   repository.
;       Modified: 11/09/11, TPEB -- Now the velocity structure element
;                                   .sig_line is calculated using the
;                                   rms_dg_spec() routine.
;       Modified: 11/28/11, TPEB -- Added error checking for ensuring
;                                   match from WHERE() statement.
;       Modified: 10/30/12, TPEB -- Migrated HCO+ / N2H+ spectra over
;                                   to Yancy's (2012) published
;                                   catalog.
;
;-

;; Funtion to return the index j_bgps which is the BGPS source for a
;;   given molecular line spectrum
FUNCTION IMPORT__J_BGPS, cnum, str, fn, nind
  
  start   = STRPOS(fn,'/',/REVERSE_SEARCH)+1
  len     = STRPOS(fn,'_',/REVERSE_SEARCH) - start
  objname = STRMID(fn,start,len)
  
  j_str = WHERE(str.name EQ objname)
  j     = WHERE(cnum EQ str[j_str].cnum, nind)
  
  RETURN, j
END

;; Function to take the input spectrum and process it into an output
;;   spectrum via interpolation and smoothing.
FUNCTION IMPORT__MAKE_SPECTRUM, inspec, hd, vstd, native_rms
  
  ;; Baseline the spectrum
  spec = baseline_spec(inspec)
  ;; EXTAST doesn't work on these headers... use homebrew code...
  v = v_from_hdr(hd)
  native_rms = rms_dg_spec(v,spec)
  
  ;; Interpolate to v_std
  specout = interpol(spec, v, vstd)
  outside = WHERE(vstd LT min(v) OR vstd GT max(v), nout)
  IF nout GT 0 THEN $
     specout[outside] = 0.
  
  ;; Use a Savitzky-Golay filter, degree 6, width = 69 points
  ;; The smoothing box is +- 6.8 km/s about the central point.
  filter = savgol(34,34,0,6)
  spectrum = convol(specout, filter, /edge_trun, /nan)
  
  RETURN, spectrum
END


;;========================================================================
;; Main Routine
PRO IMPORT_DENSEGAS_SPECTRA, PLOT=plot, HCOP=hcop, N2HP=n2hp, NH3=nh3, $
                             CS21=cs21, ORDER=order, CREATE=create
  
  COMPILE_OPT OBSOLETE
  
  ;; Set flags for plotting particular spectral sets
  plot = KEYWORD_SET(plot)
  hcop = KEYWORD_SET(hcop)
  n2hp = KEYWORD_SET(n2hp)
  nh3  = KEYWORD_SET(nh3)
  cs21 = KEYWORD_SET(cs21)
  
  ;; Set polynomial baselining order
  IF n_elements(order) EQ 0 THEN order=15
  
  ;; Check to see if velocity structure already exists.
  vstr = './local/bgps_velocity_struct.sav'
  IF ~FILE_TEST(vstr) OR KEYWORD_SET(create) THEN create_velocity_structure
  
  ;; Get current parameters
  defsysv, '!MW', exists = exists
  IF NOT exists THEN galactic_params 
  
  ;; Read in BGPS Distance Database
  s = read_bgps_csv('bgps_distance_database.csv',csv,/ver)
  
  ;; Restore all the necessary IDL structures containing
  ;;    identification information
  restore,vstr,/ver
  restore,'./local/BGPS_molec.sav',/ver
  restore,'./local/gbt_vel_fits.sav',/ver
  restore,'./local/ARO_ASSOC.sav',/ver
  restore,'./local/bgps_grs_spectra.sav',/ver
  
  ;; Define directory locations for the various dense gas tracers
  ARO_DIR = './local/aro12m/'
  
  ;; Create lists for each dense gas tracer
  ind_hcop = WHERE(molec.hcop.flag GT 0, n_hcop)
  ind_n2hp = WHERE(molec.n2hp.flag GT 0, n_n2hp)
  hht_h = molec[ind_hcop]
  hht_n = molec[ind_n2hp]
  
  n_nh3 = n_elements(gbt)
  cs21_list = FILE_SEARCH(ARO_DIR+'G*.fits', COUNT=n_cs21)
  
  print,n_hcop,n_n2hp,n_nh3,n_cs21

  velocity.n_hcop = 0
  velocity.n_n2hp = 0
  velocity.n_cs21 = 0
  velocity.n_nh3  = 0
  
  ;;=======================================================================
  ;; First, drop in the GRS 13CO spectra from the GRS structure
  message,'Copying over GRS 13CO spectra...',/inf
  velocity.grs13co = grs.spectrum
  
  
  ;;=======================================================================
  ;; Loop through HCO+ entries in Yancy's (2012) table,
  ;;   inserting gaussian fits into the velocity strucutre.
  message,'Now doing HCO+...',/inf
  
  FOR i=0L, n_hcop-1 DO BEGIN
     
     ;; Get BGPS index for this spectrum
     j_bgps = WHERE(s.cnum EQ hht_h[i].cnum, nind)
     
     ;; Check that cnum NE 0
     IF nind NE 0 THEN BEGIN
        v_std = velocity[j_bgps].v_std
        
        spectrum = gauss_1(v_std, $
                          [hht_h[i].hcop.tmb,hht_h[i].hcop.vlsr,hht_h[i].hcop.lw])
        velocity[j_bgps].sig_hcop = hht_h[i].hcop.sigt
        
        ;; Do plotting, if specified
        IF plot AND hcop THEN BEGIN
           plot,v_std,spectrum,xr=[-100,200],yr=[-0.5,4.0],/yst,$
                tit='HCO!u+!n',thick=2
           wait,0.2
        ENDIF
        
        ;; Place this information into the IDL structure
        velocity[j_bgps].n_hcop += 1
        velocity[j_bgps].hcop = spectrum
        velocity[j_bgps].has_rv_types = $
           set_rv_flag(velocity[j_bgps].has_rv_types,'hcop')
        
     ENDIF
  ENDFOR
  
  
  ;;=======================================================================
  ;; Loop through N2H+ spectra, inserting them into the velocity structure
  message,'Now doing N2H+...',/inf

  FOR i=0L, n_n2hp-1 DO BEGIN
     
     ;; Get BGPS index for this spectrum
     j_bgps = WHERE(s.cnum EQ hht_n[i].cnum, nind)
     
     ;; Check that cnum NE 0
     IF nind NE 0 THEN BEGIN
        v_std = velocity[j_bgps].v_std
        
        spectrum = gauss_1(v_std,[hht_n[i].n2hp.tmb,$
                                 hht_n[i].n2hp.vlsr,hht_n[i].n2hp.lw])
        velocity[j_bgps].sig_n2hp = hht_n[i].n2hp.sigt
        
        ;; Do plotting, if specified
        IF plot AND n2hp THEN BEGIN
           plot,v_std,spectrum,xr=[-100,200],yr=[-0.5,4.0],/yst,$
                tit='N!d2!nH!u+!n',thick=2
           wait,0.2
        ENDIF
        
        ;; Place this information into the IDL structure
        velocity[j_bgps].n_n2hp += 1
        velocity[j_bgps].n2hp = spectrum
        velocity[j_bgps].has_rv_types = $
           set_rv_flag(velocity[j_bgps].has_rv_types,'n2hp')
        
     ENDIF
  ENDFOR
  
  
  ;;=======================================================================
  ;; Loop through NH3 spectra, inserting them into the velocity
  ;; structure... these are the Gaussian fits from Miranda used in her
  ;; 2 NH3 papers (inner & outer galaxy).
  
  message,'Now doing NH3...',/inf
  
  FOR i=0L, n_nh3-1 DO BEGIN
     
     j_bgps = WHERE(s.cnum EQ gbt[i].cnum, nind)
     
     IF nind NE 0 THEN BEGIN
        v_std = velocity[j_bgps].v_std
        
        spectrum = gauss1(v_std, [gbt[i].vlsr, gbt[i].lw, gbt[i].tmb], /PEAK)
        velocity[j_bgps].sig_nh3 = 98.d-3   ;; From \citet{Dunham:2010c}
        
        ;; Do plotting, if specified
        IF plot AND nh3 THEN BEGIN
           plot,v_std,spectrum,xr=[-100,200],yr=[-0.5,4.0],/yst,tit='NH3',$
                thick=2
           wait,0.2
        ENDIF
        
        ;; Place this information into the IDL structure
        velocity[j_bgps].n_nh3 += 1
        velocity[j_bgps].nh3 = spectrum
        velocity[j_bgps].has_rv_types = $
           set_rv_flag(velocity[j_bgps].has_rv_types,'nh3')
        
     ENDIF
  ENDFOR
  
  
  ;;=======================================================================
  ;; Loop through CS2-1 spectra, inserting them into the velocity structure
  message,'Now doing CS(2-1)...',/inf
 
  ;; Check for duplicated spectra
  dup_aro = dup_elements(aro.cnum, COUNT=count, REPEATS=repeats)
  indz = WHERE(dup_aro NE 0, nind)
  IF nind NE 0 THEN BEGIN
     dup_aro = dup_aro[indz]
     repeats = repeats[indz]
     message,'Duplicate CS21 BGPS #'+strjoin(string(dup_aro,format="(I0)"),$
                                             ' ',/single),/inf  
  ENDIF
  
  FOR i=0L, n_cs21-1 DO BEGIN
     
     hdr     = headfits(cs21_list[i],/silent)
     objname = sxpar(hdr,'OBJECT')
     
     j_aro  = WHERE(aro.name EQ objname)
     j_bgps = WHERE(s.cnum EQ aro[j_aro].cnum, nind)
     
     
     ;; Check that cnum NE 0
     IF nind NE 0 THEN BEGIN
        
        v_std = velocity[j_bgps].v_std
        undefine,spec
        spec = readfits(cs21_list[i],hdr,/silent)
        
        ;; Turn input spec into output spectrum via interpolation and
        ;;   filtering 
        spectrum = import__make_spectrum(spec, hdr, v_std, rms)

        ;; Do plotting, if specified
        IF plot AND cs21 THEN BEGIN
           plot,v_std,spectrum,xr=[-100,200],yr=[-0.5,4.0],/yst,tit='CS21'
           oplot,v,spec,color='c0ff00'x
           oplot,v_std,spectrum - specout + 3., color='a0a0ff'x
           wait,0.2
        ENDIF
        
        ;; Place this information into the IDL structure
        ;; Check if this source has a duplicated spectrum...
        velocity[j_bgps].n_cs21 += 1
        indch = WHERE( dup_aro EQ s[j_bgps].cnum, n_ch )
        IF n_ch NE 0 THEN BEGIN
           repnum = ( double(repeats[indch]) )[0]
           velocity[j_bgps].cs21     += (spectrum / repnum) 
           velocity[j_bgps].sig_cs21 += (rms / repnum)
           velocity[j_bgps].cs21_fn  = velocity[j_bgps].cs21_fn+' '+cs21_list[i]
        ENDIF ELSE BEGIN
           velocity[j_bgps].cs21     = spectrum
           velocity[j_bgps].sig_cs21 = rms
           velocity[j_bgps].cs21_fn  = cs21_list[i]
        ENDELSE
        
        ;; Set the RV flag
        velocity[j_bgps].has_rv_types = $
           set_rv_flag(velocity[j_bgps].has_rv_types,'cs21')
        
     ENDIF
  ENDFOR
  
  save,velocity,filename='./local/bgps_velocity_struct.sav',/ver
  
  ;; Clean up the memory, dude
  undefine,velocity
  undefine,grs
  undefine,hht
  undefine,aro
  undefine,gbt
  
END
